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1.  SUMMARY 


Significant  progress  has  been  made  over  the  last  several  decades  towards 
the  description  and  understanding  of  atmopsheric  turbulence,  particularly  in 
the  unstable  surface  boundary  layer.  Much  of  this  progress  is  due  to  the 
development  of  simi larity^jjeory  arguments  first  postulated  by  the  Soviet 
scientist  Obukhov  in  1946.  This  report  is  an  introduction  to  basic  similarity 
theory  arguments  and  the  application  of  the  theory  to  practical  problems  in 
diffusion.  It  is  intended  to  assist  individuals  with  only  a  passing  acquain¬ 
tance  with  diffusion  theory  in  the  planning  and  execution  of  diffusion  and 
hazard  studies.  } 


fundamental  diffusion  problem  is  the  dilution  that  results  from  turbu-  • 
lent  vertical  cloud  growth.  Vertical  cloud  growth  is  described  by  the 
standard  deviation  of  the  vertical  distribution  of  cloud  material  ((o^). 
Estimates  ofuj^yare  frequently  obtained  from  Pasqui  11 -Turner  (P-T)  stability 
categories,  srgma  z  can  also  be  calculated  as  the  product  of  travel  distance 
(x)  and  the  standard  deviation  of  the  vertical  wind  angle  (fae).  An  explicit 
functional  relationship  between (oe  and  the  stability  parameter  z/L,  where  z 
is  height  and  L  is  the  characteristic  length,  is  derived  from  similarity 
theory.  This  functional  relationship  is  verified  on  a  set  of  precisely 
measured  mlcrometeorological  field  data.  The  performance  of  similarity  theory 
calculations  and  P-T  estimates  is  subsequently  tested  using  oe  data 
'rbHscted  at  DPG.  The  similarity  theory  computations  are  shown  to  be  substan¬ 
tially  superior  to  the  P-T  oe  estimates.  Adoption  of  similarity  theory 
arguments  should  materially  Improve  Army  testing  and  hazard  prediction  efforts 
involving  atmospheric  diffusion  aind  transport  processes. 


2.  INTRODUCTION 


Sigma  z  is  a  measure  of  the  vertical  dilution  in  a  diffusing  cloud  caused 
by  turbulent  motions.  Diffusion  calculations  are  critically  dependent  on  an 
accurate  specification  of  az.  In  diffusion  models  az  is  usually  speci¬ 
fied  as  the  product  of  travel  distance  (x)  and  ae,  where  ae  is  a  measure 
of  the  rate  of  vertical  mixing  (turbulence).  Sigma  e  is  defined  as  the 
vertical  wind  angle  standard  deviation.  DPG  uses  bi-directional  wind  vanes 
(bivanes)  to  measure  the  wind  angles  from  which  oe  is  calculated. 


When  measured  wind  angle  data  are  unavailable,  oz  estimates  can  be 
obtained  from  stability  categories.  A  fixed  value  of  ae  is  implied  for  each 
category.  Pasqulll  (1961)  introduced  the  most  widely  used  set  of  stability 
categories.  Turner  (1964)  modified  the  Pasquill  method  so  that  sun  elevation 
angle,  wind  speed,  and  cloud  cover  can  be  used  to  determine  the  stability 
category.  The  Pasqulll -Turner  (P-T)  method  is  popular  because  it  requires  a 
minimum  of  meteorological  data  and  a  simple  set  of  look-up  tables  to  arrive  at 
an  estimate  of  o2.  The  P-T  method  is  widely  used  for  hazard  calculations. 

It  Is  also  used  to  categorize  stability  for  smoke  diffusion  trials  and  as  a 
go/no-go  test  criterion. 


An  alternative  to  direct  ae  measurement  or  stability  category  estimates 
for  diffusion  problems  within  the  surface  boundary  layer  (SBL)  is  presented  in 
this  report.  The  SBL  Is  the  thin  layer  of  air  extending  several  tens  of 
meters  above  the  earth's  surface.  Turbulent  motions  with  the  SBL  are  deter¬ 
mined  by  fluxes,  that  is  rate  of  transfer  of  momentum  (represented  by  shear 
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stress,  t)  and  heat  (F^).  Shear  stress  and  heat  flux  are  not  directly 
measured  at  DPG  because  high-frequency  instrumentation  is  unavailable.  These 
fluxes  may  be  calculated,  however,  using  gradients  obtained  from  profile 
measurements  of  wind  speed  and  temperature.  It  Is  therefore  possible,  under 
certain  conditions,  to  specify  turbulent  motions  within  the  SBL  as  a  function 
of  appropriate  combinations  of  these  gradients. 

Similarity  theory  rises  from  the  fact  that  within  the  first  several  tens 
of  meters  above  ground  level  there  exists  a  set  of  characteristic  parameters 
with  units  of  length,  velocity,  and  temperature  which  are  essentially 
Invariant  with  height.  These  characteristic  parameters  define  the  SBL  as  a 
layer  of  nearly  constant  stress  and  heat  flux.  Monin  and  Obukhov  (1954)  used 
dimensionless  combinations  of  these  characteristic  parameters  to  describe  the 
SBL.  Of  principal  importance  in  this  report  is  the  characteristic  length  (L). 
The  dimensionless  ratio  of  height  (z)  to  L  forms  the  fundamental  stability 
property  of  the  SBL.  Consequently,  SBL  turbulence  can  be  described  as  a 
function  of  z/L,  using  similarity  theory  arguments  originated  by  Obukhov  in 
1946  (see  Obukhov,  1971),  and  sigma  e  can  be  calculated  from  measured 
gradients  of  wind  and  temperature. 

3.  PURPOSE 

The  purpose  of  this  report  is  to  derive  the  stability  parameter  z/L  as  an 
explicit  function  of  measureable  wind  and  temperature  gradients  and  to  apply 
z/L  to  the  calculation  of  ae. 

4.  SYMBOLS 

cp  specific  heat  at  constant  pressure  (calories  gram-1  degree-1) 

E  evaporation  fraction,  the  ratio  F-|/(Fi  +  Ffo)  (dimensionless) 

FE  fractional  error  (dimensionless) 

FErms  root  mean  square  fractional  error  (dimensionless) 

f  Coriolis  parameter  (second-1) 

Fg  heat  flux  into  the  ground  (milliwatts  meter-2) 

Fh,  Fi  sensible  heat  and  latent  heat  (milliwatts  meter-2) 

Fs  flux  of  substances  (variable  dimensions) 

k  von  Karman's  constant  (dimensionless) 

Kh,  Km  eddy  dlffusivity  for  heat  and  momentum  (meters2  second-1) 

K$  eddy  diffuslvlty  for  unspecified  substances  (meters2  second-1) 

L  Obukhov  length  (meters) 

l  mixing  length  (meters) 
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>  pressure  (millibars) 

R  net  radiation  (milliwatts  meter-2) 

Rlc  critical  Richardson  number  (dimensionless) 

S  a  generalized  property  in  the  air  (variable  dimensions) 

S*  a  scale  with  dimensions  of  S 

T  time-averaged  temperature  at  a  specific  height  (degrees  Kelvin) 

H  time-averaged  wind  speed  at  a  specific  height  (meter  second-1) 

u,v  wind  speed  along  the  x,  y  axis  (meters  second-1) 

z  height  of  measurement  (meters) 

z0  roughness  length  (meters) 

a  ratio  of  to  Km  (dimensionless) 

es,  ej)  rates  of  turbulent  enerqy  addition  by  shearing  and  buoyancy 
forces  (meters  second-3) 

0  potential  temperature  (degrees  Kelvin) 

v*  friction  velocity  (meters  second-1) 

p  density  (grams  centimeter-3) 

oe  vertical  wind  angle  standard  deviation  (degrees  or  radians) 

3*  square  root  of  variance  in  vertical  wind  velocity  (meters 

second-1) 

az  standard  deviation  of  vertical  distribution  of  diffusing 

material  (meters) 

t  shear  stress  (gram  centimeter-1  second-2) 

4m.  th  dimensionless  wind  shear  and  temperature  gradient 
<J>$  dimensionless  universal  function  of  stability 

<l>  diabatic  influence  function  (dimensionless) 

5.  BACKGROUND 

The  Reynolds  number  is  a  dimensionless  number  expressing  the  ratio  of 
inertial  to  viscous  forces  in  a  moving  fluid.  When  the  Reynolds  number  is 
large,  flow  is  considered  to  be  fully  turbulent.  This  equates  to  thorough 
mixing  within  the  fluid.  The  assumption  of  a  large  Reynolds  number  through 
the  SBL  permits  a  description  of  turbulent  processes  within  this  layer  using 
semi -empirical  equations  developed  from  fluid  (dynamics  laboratory  experiments. 
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Turbulent  motions  in  a  fluid  cause  the  net  movement  of  a  diffusing 
material  down  a  gradient  at  a  rate  proportional  to  the  gradient.  The 
relationship  of  the  flux  of  a  diffusing  material  to  its  gradient  is  described 
by  a  semi -empirical  gradient  transfer  equation,  equation  (1). 

Fs  =-pKs(9S/3z) 

The  gradient  transfer  equation  describes  the  flux  (Fs)  of  substance  S  as  a 
function  of  density  (p),  its  eddy  diffusivity  (Ks)  which  specifies  the 
rate  of  turbulent  mixing,  and  the  vertical  gradient  of  S  (aS/az).  For  heat, 
this  equation  becomes 


Fh  =-pKhcp (30/3z)  ^2) 

where  is  the  sensible  heat  flux,  Kh  Is  thermal  diffusivity,  Cn  Is 
specific  heat  at  constant  pressure,  and  ae/az  is  the  lapse  rate  (gradient)  of 
potential  temperature. 

Classical  mixing  length  theory  describes  the  flux  to  gradient  relationship 
for  momentum,  with  equation  (3). 


x*pt2(du/dz)|du/dz|  (3) 

Equation  (3)  presents  t  In  fluid  flowing  past  a  stationary  boundary  as  a 
function  of  velocity  gradient  and  a  mixing  length  (z)  proportional  to  distance 
from  the  boundary.  Separation  of  du/dz  terms  allows  equation  (3)  to  satisfy 
the  requirement  that  x  changes  sign  with  (du/dz).  From  equations  (1)  and  (3), 
eddy  diffusivity  for  momentum  (K^)  takes  the  form 

Km  -  A2  (du/dz)  (4) 

Eddy  and  thermal  dlfflslvlty  have  Identical  signs  and  units.  The  ratio 
Kh/Km  a  function  of  stability  and  Is  represented  by  a. 

*h/Km  »  «  (5) 

The  ratio  of  x  to  density  Is  conveniently  represented  by  friction  velocity 
(v*),  as  defined  by  equation  (6). 


v*2  ■  x/p  (6) 

Then  equation  (3)  can  be  expressed  in  terms  of  v*,  as  In  equation  (7). 

v*2  *  t2(du/dz)|du/dz|  *  Km | du/dz |  (7) 

The  heat  flux  equation  (equation  2)  can  be  further  simplified,  using 
equation  (5)  to  eliminate  Kh  and  equation  (7)  to  eliminate  i2  In  favor  of 
v*2.  The  result  Is  equation  (8),  which  describes  heat  flux  as  a  function  of 
measureable  potential  temperature  and  wind  speed  gradients. 

Fh  =  -apcpv*2(ae/az|du/dz|  (8) 
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V  w. 


It  Is  necessary  to  measure  wind  and  temperature  over  identical  time  and 
height  intervals  for  the  gradients  used  in  equation  (8)  and  all  subsequent 
gradient  equations. 

Based  on  dimensional  arguments,  Monin  and  Obukhov  (1954)  introduced  the 
general  relationship 

dS/dz  =  S*<f>s/kz  (9 

where  S  is  the  mean  amount  of  a  property  in  a  unit  mass  of  air  and  S*  is  a 
characteristic  scale  with  dimensions  of  S  and  defined  by  Fs/pv*.  The 
variable  <|>s  is  a  dimensionless  universal  function  of  stability  which  assumes 
a  value  of  unity  for  neutral  conditions.  The  von  Karman  constant  (k)  is  a 
constant  of  proportionality  (see  Businger,  1973).  For  wind  speed  gradient, 
equation  (9)  becomes 

du/dz  =  v*<j>n,/k(z+z0)  (10 

Roughness  length  (z0)  is  introduced  into  equation  (10)  to  prevent  infinite 
shear  at  the  surface.  This  equation  is  therefore  used  only  for  z>z0.  The 
term  ^  is  a  dimensionless  wind  shear  [see  equation  (21)  for  an  empirical 
expression  of  <j>m]. 


(11). 


Integration  of  equation  (10)  yields  the  logarithmic  wind  profile,  equation 


u/v*  =  (ln((z+z0)/z0)-i|>)/k 


The  magnitude  of  z0,  <0.1  m  for  open  fields,  is  sufficiently  small  enough  that 
it  is  usually  ignored  In  the  numerator  of  equation  (11).  The  diabatic  influ¬ 
ence  function  (<|>)  accounts  for  the  effects  of  thermal  stratification  on  the 
wind  profile.  Equation  (12)  presents  the  unstable  case  diabatic  influence 
function  was  derived  as  a  unique  function  of  ^  by  Paulson  (1970). 


*  -  21n((l+*,„-l)/2)  +  ln( (l+$m"2)/2)  -  2tan- 


14>m” 1  +  */2 


Obukhov  (1971)  reasoned  that  the  mean  turbulence  in  the  SBL  is  determined 
by  only  four  parameters:  density,  a  buoyancy  parameter  (g/T),  turbulent  shear 
stress,  and  vertical  heat  flux.  These  four  parameters  combine  to  form  the 
Obukhov  length  (L),  equation  (13). 

-L  =  cpPTv*3/kgFh  =  cppTT3/2Agp1/2Fh  (13) 

The  sign  of  L  is  chosen  so  that  L  is  negative  for  unstable  thermal  stratifica¬ 
tion.  Hereafter  In  this  report,  the  Obukhov  length  proceeded  by  a  minus  sign 
(-L)  Implies  heat  flux  directed  away  from  the  surface,  and  a  lapse  of  poten¬ 
tial  temperature.  With  Cp,  t,  k,  g,  and  Fh  remaining  more  or  less  con¬ 
stant, -L  will  decrease  with  height  as  T/p1/2.  Because  the  lapse  of  T/p1/2 
over  a  few  tens  of  meters  of  height  is  a  small  fraction  of  the  absolute  sur¬ 
face  value,  -L  also  remains  nearly  constant  within  the  SBL.  The  Obukhov 
length  has  units  of  meters.  Therefore  z/-L  is  dimensionless.  Obukhov  intend¬ 
ed  for  a  ratio  of  z/-L  equal  to  unity  to  Imply  a  balance  between  buoyant  ener¬ 
gy  production  and  shear  production,  and  defined  L  as  the  height  of  the  sub¬ 
layer  of  dynamic  turbulence.  It  is  shown  below  that  this  balance  is  actually 
achieved  at  a  considerably  smaller  value  of  z/-L. 
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Because  it  depicts  the  relative  contributions  of  the  buoyancy  term,  heat 

flux,  and  shear  stress  to  turbulence,  z/-L  functions  as  a  dimensionless 

stability  parameter.  When  buoyant  forces  are  negligible  and  is  small, 
the  magnitude  of  -L  is  large  with  respect  to  any  height  within  the  SBL,  and 
z/-L  is  close  to  zero.  Turbulence  for  this  condition  is  mechanical,  or  gen¬ 
erated  by  air  blowing  over  roughness  elements.  As  the  ground  is  heated,  F^ 
assumes  a  nonnegligible  valve  and  buoyant  forces  make  a  contribution  to  turbu¬ 
lence.  Monin  and  Yaglom  (1971)  found  buoyancy  effects  contributing  to  turbu¬ 
lence  as  z/-L  exceeded  0.03.  With  increasing  F^  (or  increasing  z  with 

constant  F^),  buoyancy  increases  until  the  point  is  reached  where  turbulence 

production  due  to  buoyancy  and  shear  forces  is  in  balance.  This  occurs  as 
z/-L  approaches  0.36.  When  z/-L  increases  beyond  0.36,  buoyancy  becomes 
progressively  more  dominant  until  a  state  of  windless  convection  is  reached. 
Windless  convection,  characterized  by  intense  surface  heating  with  light  and 
variable  winds,  is  the  condition  where  shear  production  and  mechanical  turbu¬ 
lence  effects  are  negligible.  This  condition  appears  as  z/-L  approaches  1.0. 

6.  DISCUSSION 

Before  specifying  the  flux-gradient  relationship  to  turbulence  with  the 
use  of  z/L,  some  underlying  assumptions  must  be  investigated.  Both  similarity 
theory  and  P-T  oe  estimates  depend  on  the  magnitude  of  surface  roughness. 
Roughness  effects  are  specified  by  the  ratio  of  height  to  roughness  length, 
ln(z/z0)  in  equation  (11).  For  an  open  field  with  roughness  of  a  few  centi¬ 
meters,  a  factor  of  2  change  in  roughness  causes  a  variation  in  u/v*  of  10 
to  15  percent  and  a  comparable  change  in  turbulence.  Therefore,  small-scale 
changes  in  roughness  have  only  a  minor  effect  on  diffusion.  This  is  fortunate 
because  roughness  is  often  difficult  to  specify  precisely.  However,  an  abrupt 
change  in  roughness,  due  to  placement  of  obstacles  in  the  flow,  can  have  a 
significant  effect  on  turbulence.  At  a  given  site,  roughness  may  vary  with 
wind  azimuth  angle,  due  to  wind  fetch  over  roughness  elements  of  different 
size.  Roughness  may  also  vary  with  wind  speed;  at  higher  wind  speeds  some 
vegetation  bends  to  reduce  drag.  When  wind  blowing  over  a  homogeneous  surface 
encounters  an  abrupt  change  in  roughness,  a  new  sublayer  of  turbulence  devel¬ 
ops.  Turbulence  characteristics  within  this  sublayer  reflect  the  influence  of 
the  new  surface  roughness.  Elliott  (1958)  described  the  height  of  this  new 
sublayer  with  the  empirical  equation 

hnew  =  °*86  x°-8z00-2  (14) 

For  equation  (14),  z0  is  the  roughness  of  the  new  surface,  x  is  the  downwind 
distance  from  the  discontinuity,  and  hnew  is  the  height  of  the  top  of  the 
new  sublayer.  For  an  obstacle  generating  a  new  roughness  of  0.1  m  placed 
100  m  upwind  of  a  meteorological  tower,  hne#  will  be  21.6  m.  Instruments  on 
the  tower  below  this  level  will  measure  turbulence  due  to  the  new  roughness 
regime  while  above  this  level  the  old  roughness  regime  predominates.  Portions 
of  a  diffusing  cloud  passing  through  different  roughness  regimes  will  diffuse 
at  different  rates,  vastly  complicating  attempts  to  analyze  or  predict  cloud 
behavior. 

The  P-T  categories  are  set  up  for  a  site  roughness  of  a  few  centimeters. 
For  maximum  accuracy,  median  turbulence  values  for  each  P-T  category  should  be 
established  using  meteorological  measurements  at  each  site  where  P-T  cate¬ 
gories  are  used.  This  was  done  at  DPG  by  Waldron  (1977).  Changes  in  rough¬ 
ness  are  accomodated  by  equation  (11)  for  similarity  theory.  For  diffusion 
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tests,  it  is  highly  desirable  to  avoid  drastic  changes  in  roughness  with  the 
resultant  formation  of  multiple  roughness  sublayers. 

The  similarity  and  P-T  methods  both  rely  on  the  assumption  of  constant 
stress  in  the  SBL  to  produce  reasonable  turbulence  estimates.  Validity  of  the 
constant  t  assumption  is  examined  using  an  equation  of  motion. 

3txz/3z  =  pdu/dt  +  3P/3x  -  pfv  (15) 

Equation  (15)  describes  the  lapse  of  t  in  the  xz  plane.  An  analogous  expres¬ 
sion  can  be  written  for  x  in  the  yz  plane.  Terms  on  the  right  side  of  equa¬ 
tion  (15)  include  the  time  rate  of  wind  speed  change  (du/dt)  in  the  alongwind 
(x)  direction,  the  x  component  of  pressure  gradient  (3P/3x),  and  the  Coriolis 
term  („fv)  where  f  is  the  Coriolis  parameter  and  v  is  the  crosswind  component 
of  wind  speed.  For  quasi -geostrophic  balance,  3P/3X  is  in  rough  equilibrium, 
with  pfv.  The  x  gradient  is  therefore  governed  largely  by  the  magnitude  of 
alongwind  acceleration  (du/dt). 

For  similarity  method  computations,  stress  is  assumed  to  be  nearly  con¬ 
stant  when  x  at  a  given  height  is  within  20  percent  of  its  surface  value,  to* 
For  v*  =  50  cm/sec  and  density  at  a  constant  value  of  0.0012  g/cm3,  equation 
(6)  yields  t0  =  3.0  g/cm  sec2.  A  strongly  accelerating  wind,  with  6  m/sec 
speed  change  within  10  minutes,  produces  a  du/dt  of  1.0  cm/sec.  From  equation 
(15),  this  acceleration  would  produce  a  Ax  of  0.6  g/cm  sec2  (20  percent  of  x0) 
at  a  height  of  5.0  m.  For  a  more  typical  wind  speed  change  of  1.0  m/sec  in  10 
minutes,  stress  remains  within  20  percent  of  x0  for  30  m.  With  extremely  low 
wind  speeds,  as  expected  in  windless  convection,  \>*  may  be  10  cm/sec  and  x0 
may  be  only  0.12  g/cm  sec2.  Then,  a  1.0  m/sec  wind  speed  change  in  10  minutes 
will  produce  a  stress  deviating  20  percent  from  x0  in  1.2  m.  Therefore,  the 
constant  x  assumption  is  valid  for  several  tens  of  meters  above  the  surface 
when  accelerations  are  modest  and  v*  is  moderately  large.  In  a  rapidly 
accelerating  wind  field,  or  a  windless  convection  regime  with  light  and 
variable  winds,  the  constant  stress  can  be  violated  within  a  few  meters  of  the 
surface.  Fortunately,  large  sustained  accelerations  in  the  wind  are  not 
likely  except  in  the  presence  of  convective  storms.  Erratic  results  can  also 
be  expected  when  light  winds  and  intense  surface  heating  persist. 

Assumptions  regarding  fluxes  of  heat  and  moisture  are  treated  differently 
by  the  P-T  and  similarity  methods.  The  P-T  method  assumes  a  constani  func¬ 
tional  relationship  between  net  solar  radiation  (R)  and  F^. 

Net  radiation,  F^,  latent  heat  flux  (F^),  and  flux  of  heat  into  the 
ground  (Fg)  are  related  by  equation  (16). 

R-Fg  =  Ft  +  Fh  (16) 

Swinbank  (1964)  tabulated  R  and  Fg  data  taken  at  an  open  field  near  Kerang, 
Australia  under  a  variety  of  weather  conditions.  These  data  show  Fg  to  be 
approximately  20  percent  of  R,  with  the  ratio  holding  fairly  constant  at  that 
site.  Changes  in  the  Fg/R  ratio  could  be  expected  for  different  ground 
cover  or  soil  types,  but  no  comprehensive  evaluation  of  this  phenomenon  is 
available. 

Dyer  (1961)  defined  evaporation  fraction  (E)  as  the  ratio  of  F]  to 
Fl+Fh,  as  in  equation  '17) 


E  =  F1/(Fi  +  Fh)  (17) 

Dyer's  data  show  that  E  varies  from  0.7  following  a  rain  event  to  0.3  on  a  dry 
day.  Consequently,  variations  from  the  assumed  relationship  between  R  and 
Fj,  can  be  as  large  as  ±40  percent. 

Because  similarity  relationships  do  not  rely  on  constancy  of  the  ratio 
Fh/(R-Fg) ,  the  effects  of  moisture  variability  are  considerably  less 
severe.  Dyer  and  Hicks  (1970)  accounted  for  latent  heat  effects  on  L  by 
modifying  equation  (13)  into  equation  (18). 

-L  =  cppTv*3/kg(Fh+0.07Fi)  (18) 

Equation  (18)  shows  that,  for  the  maximum  value  of  F]  where  the  ratio  of 
F^/Fi  is  3/7,  the  contribution  of  F]  is  only  16  percent  of  F^.  If 
available,  the  F]  factor  should  be  considered  in  length  scale  computation, 
but  omission  would  generally  result  in  an  error  of  16  percent  or  less  (except 
over  very  wet  surfaces). 

Another  implied  assumption  required  to  obtain  reasonable  estimates  of  ae 
is  that  fluxes  remain  relatively  invariant  with  respect  to  time.  Profiles  of 
temperature  and  wind  speed  adjust  slowly  to  changes  in  the  fluxes  of  these 
parameters.  Dyer  (1963)  calculated  the  time  required  for  profiles  to  reach  90 
percent  adjustment  to  the  new  condition  following  an  abrupt  change  in  surface 
flux.  The  time  required  was  14  minutes  for  profiles  to  stabilize  with  the  new 
condition  to  a  height  of  10  m,  and  30  minutes  to  a  height  of  20  m.  He  con¬ 
cluded  that  careful  selection  of  meteorological  conditions,  particularly  with 
respect  to  insolation  and  wind  fields,  is  necessary  to  satisfy  the  steady 
state  assumption  implicit  in  all  flux-gradient  considerations. 

7.  METHOD 

When  site  and  meteorological  conditions  satisfy  the  assumptions  discussed 
in  the  previous  section,  flux-gradient  relationships  can  be  used  to  calculate 
z/L.  The  flux  form  of  the  stability  parameter  is  expressed  by  equation  (19). 

z/L  *  -kgFhz/pepTv*3  (19) 

Equation  (19)  can  be  transformed  into  profile  form  for  the  unstable  case  by 
introducing  equation  (8)  for  F^  and  equation  (11)  for  v*.  After  simplifi¬ 
cation,  the  result  is  equation  (20). 

z/-L  =  -agz(39/3z/|du/dz| )[ln(z/z0)-^]/Tu  (20) 

Equation  (20)  has  several  advantages  as  a  means  of  describing  stability  in 
the  SBL.  First,  equation  (20)  does  not  require  specification  of  k.  The  abso¬ 
lute  value  of  von  Karman's  constant  is  not  precisely  known  (see  Businger, 
1973);  accepted  values  range  from  0.34  to  0.41.  Second,  equation  (20)  does 
not  contain  higher  order  terms,  such  as  (du/dz)2  as  found  in  some  other  stabi¬ 
lity  parameters.  Gradients  of  wind  and  temperature  are  small,  often  the  same 
order  of  magnitude  as  the  least  significant  digit  in  the  profiles  from  which 
they  are  computed.  Using  gradients  as  higher  order  terms  magnifies  the  poten¬ 
tial  error.  Accuracy  of  the  gradient  measurements  required  in  equation  (20) 
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represent  the  largest  source  of  error  in  the  expression.  This  is  discussed 
further  in  Verification  (Part  8)  below. 

Equation  (20)  contains  two  tehms,  a  and  41,  which  “equire  further 
specification  as  functions  of  z/-L.  The  term  a  (the  ratio  of  thermal  to  eddy 
diffusivity)  is  identified  by  Businger  et.  al .  (1971)  as  the  ratio  4>m/4>h- 
Businger,  et  ^1_.  established  empirical  relationships  between  <j>m,  <j>b»  and 
z/-L  as  presented  in  equations  (21)  and  (22). 

*m  «  (1  +  15Z/-L)-1/4  (21) 

4>h  =  0.74(1  +  9Z/-L)"1/2  (22) 

Then  u  can  be  expressed  as  a  function  of  z/-L,  as  in  equation  (23). 

a  =  1.35(1  +  9z/-L)1/2/(l  +  15z/-L) V4  (23) 

With  <t>m  described  above  as  a  function  of  z/-L,  41,  given  by  equation  (12)  as 
a  function  of  is  also  a  function  of  z/-L.  With  a  and  41  specified  as 
functions  of  z/-L,  a  simple  convergence  scheme  allows  for  rapid  computer 
solutions  of  z/-L  from  equation  (20)  using  wind  and  temperature  gradient  data 
(Appendix  B)  . 

With  z/-L  determined  from  measured  wind  and  temperature  profiles,  it  is 
possible  to  move  towards  a  similarity  theory  solution  for  ae.  Sigma  e  is 
approximated  by  the  ratio  of  ow  (the  square  root  of  the  variance  in  vertical 
wind  velocity)  and  u  as  in  equation  (24). 


<re  *  ty/u  (24) 

Panofsky  and  McCormick  (1960)  assumed  that  ow  is  determined  by  z  and  the 
rates  of  turbulent  energy  addition  by  shearing  (es)  and  buoyant  (eb) 
forces,  as  in  equation  (25). 

ow  =  A[z(es  +  Beb)]1/3  (25) 

A  and  B  in  equation  (25)  are  constants  to  be  determined  below.  Pasquill 
(1974)  defines  es  and  eb  as  follows: 

es  =  v*2d<J/dz  (26) 

eb  =  gFh/PcpT  (27) 

Equation  (13)  can  be  used  to  simplify  equation  (27),  producing  equation  (28). 

eb  =  -v*3/Lk  (28) 

Equation  (10)  is  then  used  to  eliminate  v*  from  equations  (26)  and  (28). 
Equation  (26)  becomes  equation  (29),  and  the  ratio  es/eb  is  given  by 
equation  (30). 


es  =  k2z2(du/dz)3/ 4>m2 


(29) 


es/ eb  =  "L<fcm/z 


(30) 
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Equations  (29)  and  (30)  are  then  applied  to  equation  (25).  After  some 
algebraic  manipulation  and  division  by  v*,  the  result  is  equation  (31). 

ow/v*  =  Ak  “ 1  / 3  ( 4>m  +  B(zZ-L))1/3  (31 

Following  suggestions  in  Pasquill  (1974),  constants  Ak-1/3  and  B  were 
assigned  values  of  1.3  and  1.73  respectively. 

The  dimensionless  wind  shear  term  in  equation  (31)  is  a  slowly  decreasing 
function  of  z/-L,  whereas  the  buoyancy  term  (Bz/-L)  is  an  increasing  function 
of  z/-L.  The  combined  effect  of  these  terms  is  that  aw/v*  remains  nearly 
constant  until  z/-L  reaches  a  value  of  0.36.  This  is  the  equilibrium  point 
where  <{,m  equals  1.73  z/-L  and  shear  and  buoyancy  forces  are  in  balance. 

Beyond  this  point,  the  buoyancy  effect  becomes  dominant  and  aw/v*  increas¬ 
es  rapidly  with  increasing  z/-L.  Using  equation  (21)  for  the  working 
expression  for  aw/v*  is  equation  (32). 

aw/v*  =  1.3(  (1  +  15z/-L)-1/4+  1.73Z/-L)1/3  (32 

Equation  (32),  used  in  a  ratio  with  the  logarithmic  wind  profile  [equation 
(11)],  satisfies  equation  (24)  and  produces  an  estimate  of  oe. 

<*/».  1.3((1  ♦  15Z/-L)1/'**  1.73Z/-L)1/3  .  _  _ 


u/v*  [ln(z/z0)  -  Ipj/k 

Equation  (33)  describes  oe  as  a  function  of  stability  parameter  z/-L, 
height,  and  roughness.  Following  Businger  (1973),  k  is  assigned  the  value 
0.35. 

8.  VERIFICATION 

a.  Kansas  1968  Field  Program.  The  similarity  method  for  ae  computation 
developed  in  previous  sections  was  written  into  a  short  computer  program 
(Appendix  B)  and  tested  on  Project  Kansas  field  data  (Izumi,  1971).  This 
project  was  conducted  in  the  summer  of  1968  at  a  field  site  in  southwest 
Kansas,  The  data  include  drag  plate  measurements  of  v*  and  research  grade 
wind  and  temperature  profile  measurements  on  a  32-m  tower.  Turbulence  data  at 
5.66,  11.31,  and  22.63  m  were  obtained  using  sonic  anemometers.  High  fre¬ 
quency  turbulent  velocity  fluctuations  and  air  temperature  fluctuations  were 
measured  20  times  per  second  at  three  levels  concurrently  with  temperature  and 
wind  speed  profile  measurements.  The  basic  sample  averaging  period  was  15 
minutes  and  sampling  rate  for  slow  response  profile  data  was  one  per  second. 
Equation  (24)  was  used  to  obtain  oe  for  each  unstable  case  averaging  period. 
Site  roughness  was  0.024  m.  Gradients  of  wind  speed  and  temperature  were 
obtained  from  measured  profiles  by  differentiating  second  order  polynomials 
In  (lnz)  fitted  to  measured  values  (Izumi,  1971). 

Verification  required  comparison  of  "measured"  ae  obtained  from  the  high 
frequency  turbulence  measurements  with  ae  computed  from  gradients  of  the 
tower  profile  data  via  equations  (20)  and  (33).  Fractional  error  (FE),  defin¬ 
ed  by  equation  (34),  was  used  for  verification. 


FE  =  (oecomp  -  oemeas)/1/2(oecomp  +  agmeas) 


Fractional  error  is  the  difference  between  observed  and  predicted  oe,  normal¬ 
ized  by  the  average  of  observed  plus  predicted  ae.  FE  computed  in  this  way 
provides  a  logarithmically  unbiased  means  of  presenting  the  fractional  differ 
ence  between  measured  and  computed  parameter  values.  Mean  FE  is  defined  by 
equation  (35). 

FE  =  £FE/N  (35 

Fractional  error  is  positive  (negative)  if  the  computed  ae  is  greater  (less) 
than  the  mesured  value.  Root  mean  square  fractional  error  (FErms)  is  defined 
by  equation  (36). 

FErms  =  (^FE2/N)1/2  (35 

Results  of  the  unstable  case  Kansas  field  data  verification  are  in  Table 


Table  1.  Mean  and  Root  Mean  Square  Fractional  Error  for  Kansas  1969  Field 
Program  Data  (Reference  1)  at  Indicated  Heights  Above  Ground  Level 


Hei ght 

(m) 


Mean 

Measured 

oe 


Mean 

Calculated 

°e 


FErms 


Number 

of 

Observations 


There  is  a  slight  tendency  for  underprediction  of  ae,  as  demonstrated  by 
the  negative  values  of  mean  FE  in  Table  1.  Fractional  errors  were  small  for 
almost  every  run.  At  5.66  m  there  were  no  cases  of  FE  >0.2.  At  11.31  m 
there  were  two  runs  with  FE  >0.2  (2.7  percent),  and  at  22.63  m  there  were  11 
runs  with  FE  >0.2  (16  percent).  In  no  run  did  the  error  approach  a  factor 
of  2,  or  FE  =  0.67. 

The  FErms  was  less  than  .10  at  the  5.66-  and  11.31-m  levels,  but  increased 
to  nearly  .15  at  22.63  m.  This  increase  in  error  at  22.63  m  was  largely  due  to 
resolution  of  the  gradients  to  only  two  significant  digits  beyond  the  decimal 
point.  At  22.63  m,  the  gradients  used  for  ae  computations  were  the  same 
order  of  magnitude  as  the  least  significant  digit.  Because  of  this  gradient 
problem,  -L  was  again  computed  using  the  5.66  m  gradients.  This  -L  was  assumed 
to  remain  constant  through  the  surface  boundary  layer,  and  turbulence  at  22.63 
m  was  then  re-computed  using  this  constant  -L.  The  result  was  a  decrease  in 
FErms  from  0.149  to  0.125.  In  spite  of  the  resolution  difficulties,  these  data 
demonstrate  the  effectiveness  of  the  similarity  method  turbulence  computations 
with  the  SBL. 


b.  PPG  Tower.  Measurements  obtained  with  research  grade  instruments, 
during  carefully chosen  meteorological  conditions  at  a  uniform  site,  demon¬ 
strated  the  potential  accuracy  for  similarity  method  turbulence  computations. 


However,  to  be  truly  useful  for  diffusion  calculations,  a  method  must  work 
under  a  wide  variety  of  locations  and  meteorological  conditions.  For  this 
reason,  the  effectiveness  of  the  similarity  and  P-T  ae  computations  methods 
were  compared  using  turbulence  measurements  from  a  meteorological  tower 
where  data  were  collected  hourly  over  an  eight-week  period.  These  data 
consisted  of  turbulence,  wind,  and  temperature  profiles  in  10-minute  averages 
for  every  data  collection  hour  at  the  DPG  West  Vertical  Grid  (Waldron,  1977). 
The  sampling  rate  was  one  per  second.  Gradients  of  wind  and  temperature  were 
taken  from  least  squares  fitted  data.  The  turbulence  measurements  were 
obtained  from  bivanes.  Kaimal  et  al .  (1964)  performed  a  comparison  of  bivane 
and  sonic  anemometer  vertical  turbuTence  component  measurements.  The  bivane 
data  were  found  to  correspond  well  with  the  sonic  data  for  the  lower 
frequencies  of  turbulence  expected  during  unstable  lapse  rates.  Weather 
observations  were  also  included  with  DPG  tower  data  to  allow  computation  of 
P-T  categories. 

Waldron  (1968,  unpublished  technical  report)  used  the  West  Vertical  Grid 
tower  data  to  relate  ae  values  to  each  P-T  category.  The  median  oe  for 
each  category  was  obtained  by  dividing  the  recommended  oz  at  100  m  by 
distance.  The  resultant  ae  in  radians  was  then  converted  to  degrees  (Table 
2)  for  use  in  the  evaluation. 


Table  2.  Median  8.0-m  Level  ae  Values  for  P-T  Stability 

Categories  from  DPG  West  Vertical  Grid  Turbulence 
Measurements. 


Stability  Category  ABC 

D 

E 

F 

Median  ae  15.5  10.08  6.7 

4.5 

2.98 

2.0 

Before  evaluation  of  the  P-T  and  similarity  methods,  the  data  In  Waldron 
(1977)  were  quality  checked.  Systematic  errors  In  the  temperature  profiles 
were  present  for  approximately  the  last  half  of  the  data.  This  was  probably 
caused  by  one  or  more  of  the  temperature  sensors  drifting  off  calibration. 

Data  containing  systematic  errors  were  eliminated,  as  were  incomplete  data 
records,  leaving  139  unstable  cases  for  evaluation.  For  tabulation,  the  data 
for  each  of  the  cases  were  grouped  by  the  measured  c*  value.  The  results  of 
verification  at  the  8.0-m  level  are  presented  In  Tables  3  and  4. 

The  DPG  tower  results  show  fractional  error,  increasing  with  measured  oe 
for  both  similarity  and  P-T  methods.  A  consistent  pattern  of  underprediction 
Is  also  present  with  both  methods.  The  similarity  method  was  In  error  by 
I FE I >0.2  for  37  trials  (27  percent  of  the  total),  while  the  P-T  method  posted 
an  equivalent  error  level  for  60  trials  (43  percent).  Four  similarity  method 
trial  calculations  were  In  error  by  a  factor  of  two  for  a  rate  of  3  percent, 
versus  12  trials  (9  percent)  for  the  P-T  method.  Most  of  the  large  errors  for 
both  methods  were  due  to  underprediction. 
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Table  3.  Mean  Observed  and  Calculated  DPG  West  Vertical  Grid  Turbulence 
Measurements  for  Indicated  Measured  oe  Ranges. 


Measured 

oe 

Range 

Number 

of 

Observations 

.  Mean 
Observed 
°e 

Mean 

Similarity 

®e 

Mean 

P-T 

oe 

<5.0 

44 

4.62 

4.76 

4.60 

5.1  -  10.0 

61 

6.53 

6.39 

5.73 

10.1  -  15.0 

14 

11.91 

9.77 

8.63 

>15.1 

20 

21.6 

18.15 

10.39 

Table  4. 

Mean  and  Root 
Vertical  Grid 

Mean  Square  Fractional  Error 
Turbulence  Measurements  for 

for  8.0-M  Level,  DPG  West 
Indicated  ae  Ranges. 

Measured 

Number 

°e 

of 

Si mi larity 

Si  mi  lari 

ty 

Range 

Observations 

FE 

FErms 

P-T  FE 

P-T  FErms 

<5.0 

44 

0.029 

0.080 

-0.006 

0.097 

5.1  -  10. 

0  61 

-0.040 

0.189 

-0.146 

0.277 

.0  14 

-0.134 

0.396 

-0.328 

0.385 

>15.1 

20 

-0.247 

0.419 

-0.689 

0.724 

All 

139 

-0.058 

0.242 

-0.198 

0.357 

9.  CONCLUSIONS 

DPG  West  Vertical  Grid  tower  grid  verification  reveals  that  the  P-T  turbu¬ 
lence  estimates  are  quite  accurate  for  near-neutral  values  of  ae,  but  that 
error  increases  with  the  magnitude  of  oe.  This  increase  in  error  is  due  in 
part  to  the  violation  of  the  P-T  method  assumptions  described  In  Section  6.  A 
measurement  of  mean  wind  speed  and  an  estimate  of  net  radiation  provide 
Insufficient  information  for  the  accurate  estimation  of  turbulence  under 
diabatic  conditions.  Another  source  of  inaccuracy  is  that  the  P-T  method 
offers  a  choice  of  only  six  possible  oe  values,  chosen  to  represent  centroid 
turbulence  values  within  each  category.  Turbulence  occurs  as  a  continuum  and 
varies  with  height  and  time  In  response  to  shear  and  buoyancy  forces.  Any 
attempt  to  estimate  turbulence  using  a  small  number  of  fixed  categories  must 
leave  room  for  error. 

The  similarity  theory  system  of  assumptions  and  equations,  as  derived  In 
this  report,  defines  the  relationship  between  turbulence  created  by  fluxes  of 
heat  and  momentum,  and  measured  gradients  of  temperature  and  wind.  The 
required  assumptions  are  considerably  less  restrictive  than  those  for  the  P-T 
method,  and  the  similarity  system  of  equations  presents  the  minimum  amount  of 
information  required  to  adequately  specify  the  relationship  between  measured 
profile  gradients  and  turbulence.  Because  gradients  decrease  rapidly  with 
height,  gradient  measurements  should  be  taken  at  low  levels,  preferably  at 
4.0  m. 
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Turbulent  estimates  derived  by  the  similarity  method  from  Kansas  1968 
Field  Program  profiles  (Izumi,  1971)  are  accurate,  matching  measured  turbu¬ 
lence  values  to  within  10  percent.  The  average  error  Increased  to  24  percent 
for  uncontrolled  meteorological  conditions  and  field  grade  equipment  used  at 
the  0P6  West  Vertical  Grid  tower.  Two  types  of  error  were  noted  in  these  data 
sets.  The  first  Is  a  systematic  underprediction  of  oe  by  the  similarity 
method  for  all  except  the  near  neutral  turbulence  values.  Increasing  constant 
B  In  equation  (36)  from  1.73  to  a  value  near  1.8  would  largely  eliminate  the 
systematic  error.  An  interesting  consequence  of  this  change  is  that  the 
balance  between  eb  and  es  would  occur  at  z/-L  =  0.35,  which  is  equal  to 
von  Kantian' s  constant.  From  equation  (19)  It  is  easily  seen  that  when  z/-L  = 
k,  gFhz/pcpTv*3  equals  1.0,  and  the  balance  between  buoyancy  and  t 
occurs  at  unity  as  Obukhov  (1971)  intended. 

The  remaining  errors  in  random  scatter  are  due  in  part  to  violation  of 
similarity  assumptions  and  in  part  to  lack  of  temperature  measurement  accuracy 
in  the  DPG  tower  data.  The  random  error  can  be  reduced  to  the  10-percent 
value  found  with  the  Kansas  field  data  (Izumi,  1971)  by  improving  Instrument 
quality  and  data  reduction  procedures,  and  by  carefully  choosing  test 
conditions. 

10.  LIMITATIONS 

The  similarity  theory  arguments  presented  In  this  report  are  valid  only 
for  neutral  through  unstable  thermal  stratification.  A  corresponding  set  of 
equations  has  been  offered  for  the  stable  case,  but  these  equations  apply  only 
to  stabilities  less  than  the  critical  Richardson  number  (Ric).  Oke  (1970) 
states  that  beyond  the  Ric,  the  atmosphere  is  not  fully  turbulent  and  there 
may  be  no  definite  forms  for  profiles  of  wind  and  temperature.  He  concludes 
that  a  purely  empirical  approach  Is  needed  for  estimating  fluxes  beyond  Ric. 
With  clear  skies  and  strong  outgoing  radiation,  nocturnal  conditons  at  DPG 
almost  always  exceed  Ric  and,  therefore,  do  not  lend  themselves  to  a  similar¬ 
ity  theory  solution. 

Similarity  arguments  apply  only  to  the  vertical  component  of  turbulence. 
Calder  (1966)  provided  a  theoretical  demonstration  to  support  observations 
that  the  similarity  theory  cannot  be  applied  to  variances  of  horizontal  compo¬ 
nents  of  wind  velocity.  This  veiwpolnt  Is  reinforced  by  Lumley  and  Panofsky 
(1964),  and  Panofsky  et  al .  (1979).  The  horizontal  components  of  turbulence 
can  be  reasonably  obtained*  only  by  measurement. 

Similarity  arguments  do  not  take  Into  account  the  mechanical  dampening  of 
turbulence  at  the  earth's  surface.  The  turbulence  spectrum  near  the  ground  is 
distorted  by  the  proximity  of  a  fixed  boundary  which  breaks  up  large  eddies. 
This  effect  Is  particularly  noticeable  for  low  frequency,  high  amplitude 
eddies  which  develop  during  unstable  thermal  stratification.  It  Is  therefore 
not  desirable  to  apply  the  similarity  theory  sigma  e  calculation  method  to 
heights  much  below  4.0  meters. 

Similarity  arguments  presented  In  this  report  apply  only  to  the  surface 
boundary  layer.  Considerable  work  has  been  done  to  extend  the  similarity 
relationship  through  the  depth  of  the  planetary  boundary  layer.  However,  this 
Involves  scaling  with  height  of  the  mixing  layer,  a  measurement  not  presently 
available  at  DPG. 


11.  RECOMMENDATIONS 


It  is  recommended  that  the  similarity  theory  method  described  herein  be 
used  in  the  following  applications: 

a.  Hazard  Prediction.  Hazard  prediction  programs  presently  require  the 
input  of  P-T  categories.  A  simple  modification  of  these  programs  for  similar¬ 
ity  method  computations  can  improve  prediction  accuracy  if  temperature  and 
wind  profile  data  are  made  routinely  available. 

b.  Test  Go/No-Go  Criteria.  Meteorologists  and  Project  Officers  can  use 
stability  information  as  a  quantitative  means  of  evaluating  the  suitability  of 
meteorological  conditions  for  diffusion  tests.  Tower  profile  data  are  pre¬ 
sently  collected  at  the  rate  of  one  per  second  by  data  acquisition  systems. 
Gradients  from  these  data  can  be  used  to  produce  a  real-time  update  of  the 
stability  parameter  z/-L  for  test  control.  For  tests  requiring  near  neutral 
conditions,  a  z/-L  less  than  0.03  is  required.  A  z/-L  less  than  or  equal  to 
0.36  may  be  used  as  a  go/no-go  criterion  for  testing  under  free  convection 
conditions.  Beyond  z/-L  =  0.36,  buoyancy  forces  become  large  and  the  proba¬ 
bility  of  erratic  cloud  behavior  increases.  As  z/-L  approaches  1.0,  a  condi¬ 
tion  of  light  and  variable  wind  may  occur  where  the  cloud  centroid  rises  off 
the  ground.  Go/no-go  criteria  for  diffusion  testing  can  be  based  on  these 
z/-L  ranges. 

c.  Test  Data  Qjality  Control.  Similarity  theory  equations  describe  the 
relationship  between  turbulence  and  measured  meteorological  profiles.  With 
high  quality  data  collected  at  a  site  of  uniform  roughness  during  quasi- 
stationary  meteorological  conditions,  the  computed  turbulence  values  should  be 
quite  accurate.  During  testing,  comparison  of  calculated  and  measured  turbu¬ 
lence  can  be  made.  If  a  discrepancy  in  excess  of  20  percent  is  found,  then 
either  the  equipment  Is  malfunctioning  or  meteorological  conditions  are  too 
variable  for  the  collection  of  a  coherent  data  set. 

d.  Modeling.  Because  L  is  constant  within  the  surface  boundary  layer, 
profile  data  used  to  calculate  L  at  one  representative  level  provides  suffi¬ 
cient  information  for  model  computation  of  z/L  over  a  full  range  of  height 
within  the  SBL.  Since  equation  (33)  describes  oe  as  a  function  of  only  z/z0 
and  z/-L,  an  entire  turbulence  profile  up  to  heights  of  20  or  30  meters  can  be 
developed  from  a  single  set  of  lower  level  gradient  measurements.  Conse¬ 
quently,  Improved  diffusion  modeling  within  the  SBL  may  be  achieved  without 
the  requirement  for  extensive  measurements  on  tall  towers. 
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APPENDIX  B.  PROGRAM  LISTING 


ICONTROL  USLINIT, NOSOURCE 

C  PROGRAM  SIGECOMP  COMPUTES  THE  UNSTABLE  CASE  DIMENSIONLESS 
C  STABILITY  PARAMETER  Z/-L  (ZL),  WHERE  Z  IS  HEIGHT  IN  METERS 
C  AND  L  IS  THE  MON  IN-OBUKHOV  LENGTH.  WIND  AND  TEMPERATURE 
C  GRADIENT  DATA  ARE  USED  AS  INPUT  FOR  THIS  COMPUTATION.  Z 
C  IS  THE  HEIGHT  ABOVE  GROUND  FOR  WHICH  WIND  AND  TEMPERATURE 
C  GRADIENTS  WERE  COMPUTED.  ZO  IS  ROUGHNESS  LENGTH  IN  METERS. 

C  THE  PROGRAM  USES  ZL  TO  COMPUTE  THE  VERTICAL  WIND  ANGLE  STAND 
C  ARD  DEVIATION  (SIGE)  AND  COMPARES  IT  TO  A  MEASURED  VALUE 
C  (SIGM)  BY  MEANS  OF  FRACTIONAL  ERROR  (FE).  PASQU I LL -TURNER 
C  STABILITY  CATEGORY  SIGMAS  (ICAT)  ARE  ALSO  COMPARED  TO  SIGM  AND 
C  EVALUATED  USING  FE. 

EXTERNAL  ATAN 
DISPLAY  "Z  =  \"ZO  =  " 

ACCEPT  Z.ZO 
10  CONTINUE 

PI=3. 141592654 
DISPLAY  “CAT * " 

ACCEPT  ICAT 
IF ( I  CAT . EQ . 0 .  )  STOP 
IF ( I  CAT . EQ  .  1 )  CSIGE  =  15 . 5 
I F ( I  CAT . EQ . 2  )  CSI GE  =  10 . 08 
IF( I  CAT. EQ . 3  )  CSIGE=6. 7 
IF ( ICAT.EQ .4  )  CSIGE  =4 . 5 
IF( I  CAT . EQ . 5 )  CSIGE=2.98 
IF ( I  CAT . EQ  .  6  )  CSIGE=2.0 
DISPLAY  "U= , T= , DUDZ= , DTDZ= " 

ACCEPT  U , T , DUDZ , DTOZ 
E=AL0G ( Z/ZO ) 

DUDZ  =ABS ( DUDZ  ) 

IF(DTDZ.GE.O.O)  GO  TO  300 
C  VK  IS  THE  VON  KARMAN  CONSTANT 
VK= . 35 

C  G  IS  FORCE  OF  GRAVITY,  M/SEC**2 
G=9 . 8 

T=T+273 . 16 

C=-G*Z*DTDZ/(T*U*DUDZ) 

C  GO  THROUGH  ITERATIVE  PROCEDURE  TO  FIND  Z/-L 
C  ZL  IS  Z  DIVIDED  BY  THE  UNSTABLE  CASE  M-0  LENGTH,  -L 
ZL*. 00001 

20  PHI=(l.+15.*ZL)**{-.25) 

P.SI=2.*AL0G{  (l.  +  l./PHI  )/2.)+AL0G(  (l.  +  l./PHI**2)/2.)-2.* 
1ATAN ( 1 . /PHI )+P 1/2 . 

ALPHA«(1.35*{1.+9.*ZL)**.5)/(1.+15.*ZL)**.25 
ZLNEW"C*ALPHA*(E-PSI ) 

IF(ABS(ZLNEW-ZL) .LT. 0.0005)  GO  TO  50 
ZL  * ( ZL+ZLNEW ) / 2 . 

GO  TO  20 

50  ZL  -  (ZL+ZLNEW ) /2 . 

PHI»(l.+15.*ZL)**(-.25) 

PSI=2.*AL0G( (1 .  +  1./PHI )/2. )+AL0G( ( 1 . +1 . /PHI**2  )/2 . )-2.* 


1 ATAN ( 1 . /PHI )+PI/2. 

USTAR«U*VK/(E-PSI ) 

UBYUSTAR=(E-PSI )/VK 

SIGWBYUSTAR=1.3*( (1. +15. *ZL )**{-. 25 )+l . 73*ZL ) ** . 33333 
SIGE=SIGWBYUSTAR/UBYUSTAR 
C  CONVERT  SIGE  FROM  RADIANS  TO  DEGREES. 

SIGE*SIGE*180/PI 
DISPLAY  "SIGE*", SIGE 
DISPLAY  "SIGE  MEAS*" 

ACCEPT  SIGM 

FE=(SIGE-SIGM)/( ,5*(SIGE+SIGM) ) 
CFE=(CSIGE-SIGM)/(.5*(CSIGE+SIGM)) 

FERMS* ( FE ) **2 
CFERMS*( CFE )**2 
60  WR I TE ( 6 , 100  ) 

100  FORMAT ( '  1 , IX  ,  '  U  '  , 5X , ' DTDZ 1 ,2X, 1 DUDZ ' , 2X , ' CAT ' , 2X , 1 CSIGE ' , 2X , 
l'SIGM' ,4X, 'SIGE' ,4X, 'FE' ,4X, 'FERMS' , 3X , ' CFE ' , 2X , ' CFERMS  '  ,  2X  , 

1 ' Z/-L ' ) 

WR ITE(6, 200  )U, DTDZ, DUDZ, IC AT, CSIGE, SIGM, SIGE, FE, FERMS, CFE, 

1 CF E RMS  Z L 

200  FORMAT!  ‘  ' ,F5.2,1X,F6.4,1X,F5.4,2X, 11 , 2X , F5 . 2 , IX , F7 . 3 , IX , F7 . 3 
11X,F7.3,1X,F5.3,1X,F7.3,1X,F5.3,1X,F7.3) 

WRITE  (6,500) 

500  FORMAT  ('O’) 

GO  TO  400 

300  DISPLAY  "DTDZ  GE  0  COMP  NOT  POSSIBLE" 

400  ACCEPT  IDUM 
GO  TO  10 
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